function [X] = rda(Training_Samples, Testing_Samples, Training_Labels, lambda, gamma)
[intervals, class_sizes, class_count] = split(Training_Labels);
test_rows = size(Testing_Samples, 1);
[training_rows, training_columns] = size(Training_Samples);
PooledCov  = zeros(training_columns, training_columns);
X = zeros(test_rows, class_count);
N = 0;

for class = 1:1:class_count,
    k = Training_Samples(intervals(class,1):intervals(class,2),:);
    ni = class_sizes(2, class);
    PooledCov = PooledCov + (ni - 1)*cov(k, 1);
    N = N + ni;
end

PooledCov = PooledCov / (N - class_count);
for class = 1:1:class_count,
    k = Training_Samples(intervals(class,1):intervals(class,2),:);
    ni = class_sizes(2, class);
    Si_lambda = ((1 - lambda)*(ni - 1)*cov(k, 1) + lambda*(N - class_count)*PooledCov)/((1 - lambda)*ni + lambda*N);
    Si_lambda_gamma = (1 - gamma)*Si_lambda + gamma*(trace(Si_lambda)/training_columns)*eye(training_columns);
    u = mean(k);

    inv_sigma = inv(Si_lambda_gamma);
    log_det_sigma = log(det(Si_lambda_gamma));

    for row = 1:1:test_rows,
        A = Testing_Samples(row:row,:);
        X(row, class) = (A - u) * inv_sigma * (A - u)' + log_det_sigma;
    end
end
end
